The effect of multiple paternity on genetic diversity during and after colonisation 
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In metapopulations, genetic variation of local populations is influenced by the genetic content of the founders, 
and of migrants following establishment. We analyse the effect of multiple paternity on genetic diversity using 
a model in which the highly promiscuous marine snail Littorina saxatilis expands from a mainland to colonise 
initially empty islands of an archipelago. Migrant females carry a large number of eggs fertilised by 1 — 10 
mates. We quantify the genetic diversity of the population in terms of its heterozygosity: initially during the 
transient colonisation process, and at long times when the population has reached an equilibrium state with 
migration. During colonisation, multiple paternity increases the heterozygosity by 10 — 300% in comparison 
with the case of single paternity. The equilibrium state, by contrast, is less strongly affected: multiple paternity 
gives rise to 10 — 50% higher heterozygosity compared with single paternity. Further we find that far from 
the mainland, new mutations spreading from the mainland cause bursts of high genetic diversity separated by 
long periods of low diversity. This effect is boosted by multiple paternity. We conclude that multiple paternity 
facilitates colonisation and maintenance of small populations, whether or not this is the main cause for the 
evolution of extreme promiscuity in Littorina saxatilis. 
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effects. 



I. INTRODUCTION 

When new local populations are established in a metapopu- 
lation, genetic variation within the newly founded populations 
is initially governed by the genetic content of founders. At a 
later stage, during continued input of variation through mi- 
gration, the genetic composition of migrants may potentially 
contribute new variation and hence counteract loss by drift 
and selection. In brooding and sexual species, empty sites are 
most likely colonised by single fertilised females that bring 
a brood of offspring along, while founders of virgin females, 
males, and juveniles fail to mate in an empty site |Q]]. In brood- 
ing species, female promiscuity (the propensity to mate mul- 
tiple males) may influence the genetic variation carried by the 
founders, and, if so, will have consequences on the effective 
population size of the new population. 

Females mating multiple males have broods of offspring 
sired by more than one male, unless sperm competition or 
cryptic female choice prevent this. Female promiscuity, once 
believed to be rare in nature, is observed in a number of 
animal species including mammals, amphibians, fishes and 
invertebrates JH-Q]. Genotyping offspring of species with 
promiscuous females shows that a large proportion of the fe- 
males releases offspring sired by 2 — 4 males 0]. In some 
species of fish and invertebrates, levels of multiple paternity 
are even higher, since the number of males siring a brood reg- 
ularly reaches 6 — 10 |@-l8t]. An extreme example of high 
multiple paternity is the marine snail Littorina saxatilis with 
15 — 23 males siring broods of single females [9[]. In this 
study we construct a mating model for this species and anal- 
yse how multiple paternity affects population genetic variation 
and structure in a metapopulation. We consider established 
populations in equilibrium, but also populations under estab- 
lishment (during initial colonisation of a previously empty 



habitat). 

Littorina saxatilis is strictly intertidal and most abundant 
in rocky shores in the north Atlantic, with population densi- 
ties of tens to hundreds of snails per square metre ifioll . In 
contrast to many other marine snails, L. saxatilis does not 
have pelagic eggs or larvae, and therefore dispersal over a 
few metres range is infrequent. However, snails occasionally 
migrate among islands, probably by rafting. It has been esti- 
mated that within an archipelago of small and large islands, 
3% of the small islands receive a migrant snail each year Hill . 
In many areas, L. saxatilis forms metapopulations with local 
populations inhabiting discrete localities, such as islands of an 
archipelago, rock y ou tcrops and breakwaters intermingled by 
sandy substrates IU21, [l3ll . During the retreat of the ice sheet 
12000 — 15000 years ago, L. saxatilis spread from refuge ar- 
eas both in the northern Atlantic and south of the ice-cap OJ. 
Part of this postglacial expansion comprised colonisation of 
hundreds of islands in the archipelago along the Swedish west 
coast that successively became available by isostatic uplift, a 
process that is still ongoing. In this system, populations on the 
mainland and large islands are the oldest and largest, and these 
are likely to act as the ultimate sources of genetic variation 
during colonisation of emerging islands in a stepwise manner 
(Fig.[T}V). We have re-analysed genetic data from L. saxatilis 
populations in the archipelago of west Sweden and found that 
the first principal component shows a largely linear relation- 
ship between population genetic variation and size/age of the 
islands/populations, with mainland populations at the one end, 
skerry populations (skerry sizes w 10m 2 ) at the other end, 
and island populations (island sizes ps 10 5 m 2 ) in the middle 
(Fig.QJi). This suggests a simple linear stepping-stone model 
with the mainland population acting as a source for colonisa- 
tion of islands at successively younger age, and at increasing 
distance from the mainland (Fig.[T]C). 
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In this paper, we investigate how multiple paternity in L. 
saxatilis affects spatial and temporal structure of neutral ge- 
netic diversity in a metapopulation of this species. We analyse 
the effective genetic population size resulting from the mating 
behaviour observed in earlier studies, and derive simple ap- 
proximations describing the genetic diversity of populations 
during colonisation, and in its equilibrium state with migra- 
tion. We use simulations to assess the temporal effects of mi- 
gration on genetic diversity as new mutations from the main- 
land spread to distant islands. 



II. DESCRIPTION OF THE MODEL 

We construct a stepping-stone colonisation model with 
the following basic assumptions mimicking how L. saxatilis 
colonises the post-glacial archipelago of western Sweden. (1) 
Colonisation is frequent and rapid, as rafted fertilised females 
release a few hundred offspring already in the first generation 
ifTTll . (2) Small skerries are likely to be colonised within a 
few years after emergence, and hence all newly established 
populations are limited by the small size of the habitat, result- 
ing in census sizes of ps 10 2 — 10 3 Hill . (3) Colonisation is 
likely to take place in a stepping-stone manner with smaller 
and more distantly related islands being colonised from their 
closest already colonised islands. For simplicity, we consider 
a system consisting of a mainland and of linearly arranged is- 
lands of equal carrying capacities (substantially smaller than 
that of the mainland). 



A. Stepping-stone colonisation model 

In our model, islands are linearly arranged and numbered 
from 1 to k, with k being furthest away from the mainland 
(see Fig.Q]C). We include high values of k in our model (such 
as k = 10) in order to be able to assess saturation effects. The 
mainland is labelled by 0. Generations are assumed to be non- 
overlapping. At the generation when the process of colonisa- 
tion starts, the mainland is the only populated habitat, and the 
population heterozygosity on the mainland is stationary (that 
is, the mainland population is assumed to be old). 

Within our model, an empty island becomes fully colonised 
in a single generation after the arrival of one or more founder 
females from the nearest neighbouring island. This is moti- 
vated by the very large capacity for population growth of L. 
saxatilis in a suitable habitat ||12ll . In our model the founder 
females give rise to 2N offspring in total, with equal sex ratio, 
where 2N denotes the carrying capacity of an island. Upon 
establishing a given island population, its population size re- 
mains constant over time. In our model, mating takes place 
before migration, which allows us to trace only the movement 
of females (males also migrate, but since they will not mate 
after migration, they do not contribute to the progeny on the 
island they migrated to). Individuals are equally likely to mi- 
grate to each of their closest neighbours (but the mainland and 
the last island have only a single neighbour, Fig.Q]C). On av- 
erage, M females migrate per generation from one island to 
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FIG. 1 : (A) Physical structure of the marine habitats of Littorina sax- 
atilis: mainland (red), islands (green), and skerries (blue). (B) Princi- 
pal components of allozyme population differentiation in L. saxatilis 
(data from fl2ll . the presumably selected locus Aat — 1 is excluded). 
Populations are classified as mainlands (red), islands (green), or sker- 
ries (blue). (C) Stepping-stone model of a section of the archipelago, 
with the mainland (labelled by 0) acting as a source for establishing 
the island populations (1 to k). 

a neighbouring island, except for empty islands that only re- 
ceive migrants. 

In addition to the above, we assume that the population size 
on the mainland is much larger than the population size of a 
colonised island (unlike other models which assume that all 
habitats have the same carrying capacity, see, for example, 
ITill ). This allows us to treat the mainland as the only source 
of genetic variation to the island populations. In our com- 
puter simulations (see Sections II-IV in the electronic supple- 
mentary material, ESM) we set the mainland heterozygosity 
to unity. This simplifies the simulations, since the dynamics 
of the mainland does not need to be simulated explicitly. 



B. Mating model 

In order to study the consequences of multiple paternity for 
genetic diversity, we introduce a mating model to describe dif- 
ferent levels of multiple paternity in mating systems. 

Based on known life-history traits fl2ll . we assume that the 
duration of the reproductive cycle of females is the same for 
all females. Each female carries beneath her shell juveniles 
of varying degrees of maturity, and juveniles are released at 
an approximately steady rate. We also assume that after a 
successful mating, the mated female obtains a sperm package 
able to fertilise female eggs during its persistence time. Our 
observations show that sperm can be stored and used up to 
a year after mating. The number of eggs fertilised by a sin- 
gle sperm package is denoted by N cggs , and we assume that 
this number is the same for all sperm packages that a female 
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receives during the reproductive cycle. The total number of 
sperm packages received by a female during her reproductive 
cycle is denoted by I. 

The probability p that two eggs are fertilised by the same 
sperm package is p = i _1 assuming that all sperm packages 
a female received during her reproductive cycle persist until 
the end of the reproductive cycle, that all eggs are fertilised 
after all sperm packages have been collected, and that sperm 
packages are chosen with replacement to fertilise eggs. 

The scheme presented above models the process of mating 
at an individual level. We assume that individuals belong to a 
well mixed diploid population of N m males and Nf females, 
and we take N m 1 and Nf 3> 1. In our model, a female 
encounters s < N m different males during her reproductive 
cycle. Having s < N m reflects the limited movement of snails 
during the reproductive cycle. To simplify the analysis, we 
assume that all males a female encounters are equally likely 
to be her mating partner in any of the matings she experiences. 
Moreover, we assume that all females are equally successful 
mothers. Within the model described, the effective size of a 
single local population is given by (see Section I of the ESM): 

Here, k is the probability that two offspring having the same 
mother share a father 



K = p+(l-p)-. (2) 
s 

Since K decreases as the number of available mates s of a 
female increases (keeping the value of p < 1 fixed), we take 
k" 1 as a measure of the degree of multiple paternity. 

Our model reduces to the model in |fl5tl in the case a female 
encounters all males present in the population (upon substi- 
tuting the number of matings in lfl5ll by If each female 
mates with all males in the population and the probability that 
two eggs are fertilised by the same sperm package is p = 0, 
our model reduces to random mating. 

We show in Fig. [2] that the effective population size in- 
creases as s increases. For the parameters set in Fig. [21 the 
maximum value of iV c is equal to N m + Nf, which corre- 
sponds to the effective population size under random mating. 
The increasing trend of N c saturates at s « 10 for a given 
value of p. In summary, by increasing the degree of multi- 
ple paternity, the effective population size becomes larger (as 
found in 111611 for a different mating model). 

We compared the male family sizes (the number of off- 
spring of a father involved in siring the brood of a given fe- 
male) obtained under our model to those obtained under ex- 
perimental conditions, as well as in natural populations. We 
conducted experiments such that 6 virgin females were placed 
in separate aquaria, and each was accompanied with s = 10 
males. The sire of each offspring produced during a year was 
determined by genotyping 8 microsatellite loci. 
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FIG. 2: Effective population size (for the given number of available 
mates, s, and the probability that two eggs are fertilised by the same 
sperm package, p) relative to 7V m + Nf. The grey region depicts the 
average number of fathers in four broods of L. saxatilis [3], Parame- 
ters: N m = Nf = 10 3 . 

III. RESULTS 

In Fig. [3]\ we show the histograms of male family sizes 
obtained experimentally and in Fig.[3jJ we show similar data 
collected from females in natural habitats |9|]. For both sets 
of data we use computer simulations in order to find the pa- 
rameters in our model resulting in male family sizes that are 
closest to those empirically observed (the brood sizes anal- 
ysed in computer simulations correspond to those from the 
empirical data). For data obtained under experimental condi- 
tions, we vary the probability p in the computer simulations, 
and for data collected in natural habitats, we vary both s and 
p. We use a x 2 -test to quantify the distance between the em- 
pirical and simulated data. The best fits obtained are shown 
in Fig. [3JV-B by circles, and they correspond to p = 1/15 
(Fig.|3jV), and to s = 20, p = (Fig. [3^). In summary, Fig. [3] 
suggests that our mating model describes empirical data ob- 
tained in the experimental setup well, whereas the agreement 
between the model and the empirical data taken from natural 
populations is less good. This is discussed in SectionllVl 

To address the question of how multiple paternity affects 
genetic variation and structure in our metapopulation, we 
analyse genetic diversity under the colonisation model de- 
scribed in Methods. We analyse separately two phases of pop- 
ulation dynamics on each island: initial colonisation, and the 
equilibrium state that develops once the colonisation phase 
is over. For a given island, we compute the expected het- 
erozygosity in the generation when the island is colonised 
(colonisation phase), as well as the heterozygosity in the equi- 
librium state. The corresponding analytical computations are 
described in detail in Sections II and III of the ESM. We also 
study the temporal changes of heterozygosity under our model 
by computer simulations. In the following two subsections we 
present separately the results for the colonisation phase and 
for the equilibrium state. 

A. Colonisation phase 

We compute the heterozygosity du ring colonisation ana- 
lytically using a coalescent approach IU7I1 . We represent the 
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FIG. 3: Histograms of male family sizes within broods of females. 
Bars in panel A show the empirical data obtained under experimental 
conditions for s — 10; the data correspond to six broods, two of size 
20, three of size 19, and one of size 16. Bars in panel B show results 
taken from QJ; the data correspond to four broods of sizes 79, 77, 
71, and 53. The width of the bins are chosen so that the expected 
number of counts in each bin is not smaller than 5. The probability 
assigned to each bar is proportional to the bar area. Symbols and 
error bars show the result of the best fit to the experimental data, 
together with their 95% confidence intervals: p = 1/15 in A, and 
p = 0, s = 20 in B. We simulated 10' ! independent realisations of 
the mating process. 



population-size history of the population on island i as a se- 
quence of i bottlenecks (i is the number of colonisation events 
that the population ancestral to that on island i went through 
before the island was colonised). Our analytical result is valid 
for small migration rate, M <C 1. Under this assumption, 
colonisation of empty islands occurs rarely, but when it does, 
an island typically receives a single founder female (see Sec- 
tion II of the ESM). The result is given in Eq. (14) of the ESM 
and in Fig. |4}V. We see that the colonisation-phase heterozy- 
gosity decays as distance from the mainland increases. We 
note that the results of our computer simulations (see Fig. 1 A- 
C in the ESM) agree well with the analytical results for low 
migration rates. For large migration rates (M = 0.5. i.e. 
0.5 females on average per generation) by contrast, the simu- 
lations assume somewhat higher values than the theory. The 
reason for this deviation is that for M = 0.5 it is probable that 
more than one founder female comes to an empty island to es- 
tablish the population, and, consequently, will contribute with 
more genetic variation than just one founder female. For nat- 
ural populations of L. saxatilis it has been estimated that 3% 
of empty islands receive a migrant each year flUl . This esti- 
mate is close to the lower value of M used in our simulations 
(M = 0.05). An important result shown in Fig. [4}V is that at 
any distance from the mainland, multiple paternity results in 
higher heterozygosity than single paternity. Mating two males 
(s = 2) increases the values of single-paternity heterozygosity 
by 10 — 100% for the parameters used in Fig. [4jV, and mat- 
ing ten males (s = 10) increases the values of single-paternity 
heterozygosity by 10 — 300% (Fig.[4}V). The largest increase 
is observed at the island furthest from the mainland. We also 
note that mating more than about 10 males only marginally 
increases the heterozygosity (results not shown), as found in 
the case of freely mixing populations (see Subsection lHBb . 



FIG. 4: Analytically computed heterozygosity during colonisation 
(A), and in the equilibrium state (B). The lines shown from top to 
bottom correspond to: s = 10, s = 5, s = 3, s = 2, and s = 1. 
Remaining parameters: the mainland heterozygosity is set to unity, 
scaled female migration rate M = 0.05, number of females in each 
populated island N = 100, probability that two eggs are fertilised 
by the same sperm package p — 0.1, number of islands k — 10. 



B. Equilibrium state 

We show in Section III of the ESM how to compute the 
heterozygosity in the equilibrium state at distance i from the 
mainland using recursion relations (see, for example, lfl8ll ). 
Note that this derivation does not require M to be small. 

The results are given in Section III of the ESM and in 
Fig. |4jJ. As in the colonisation phase, the equilibrium-state 
heterozygosity decreases as distance from the mainland in- 
creases. Also, by increasing the degree of multiple paternity, 
the heterozygosity at a given island increases (this effect sat- 
urates at s 10, results not shown). In contrast to the strong 
effect of multiple paternity during colonisation, the effect is 
substantially smaller in the equilibrium state. We find that 
the single-paternity heterozygosity in the equilibrium state in- 
creases by 10 - 20% for s = 2, and by 20 - 50% for s = 10 
(Fig. |4j5). As in the colonisation phase, the largest increase is 
observed at the island furthest from the mainland. 

In addition, we examined the variation in heterozygosity 
over consecutive generations in a particular realisation of our 
model. We find that the heterozygosity shows strong temporal 
fluctuations (Fig. [5}V). Notably, the fluctuations are strongest 
furthest from the mainland, with periods of high diversity sep- 
arated by long periods of near or complete fixation. Hence the 
distribution of heterozygosity at large distance from the main- 
land is bimodal. The heterozygosity is expected to have a 
bimodal distribution in the case of a very small rate of income 
of new genetic material, as pointed out in lfl9ll (see Fig. 1 in 

My 

In what follows, we analyse how the durations of the phases 
of low and high heterozygosity are affected by multiple pater- 
nity. Using the results of computer simulations, we compute 
the average durations of low- and high-heterozygosity phases 
at the island furthest from the mainland. We also derive cor- 
responding analytical results under the assumption that the 
scaled migration rate M is small (see Section IV of the ESM). 
For island i = 10 we show in Fig. |5jJ the durations of low- 
and high-heterozygosity phases relative to their corresponding 



5 




0.9-1 — , — , — i i i i i— i — i 0.9^ — , — , — , — , — , — , — , — , — i 
123456789 10 123456789 10 



Number of mates, s Number of mates, s 



FIG. 5: (A) Heterozygosity as a function of the distance from the 
mainland and of time (single realisation of the model described). 
Mainland is not shown. The data correspond to 10 J generations af- 
ter the initial 7 ■ 10 6 generations. The number of available mates 
is s = 10. (B) Durations of low-and high-heterozygosity phases 
(blue, and red) relative to their corresponding values for s = 1. C 
Equilibrium-state heterozygosity (black) relative to its corresponding 
value for s = 1. Remaining parameters are as in Fig.|4j5. 



values for a single mate (s = 1). Fig. [5ji shows that multi- 
ple paternity prolongs the duration of the high-heterozygosity 
phase, and decreases the duration of the low-heterozygosity 
phase. For example, the high-heterozygosity phase for the 
highest level of multiple paternity shown (s = 10) is pro- 
longed by around 40% compared to its value under single pa- 
ternity (s = 1). The low-heterozygosity phase is shortened 
by only around 10% for s = 10 (Fig. |5j$). For comparison, 
Fig-SC shows the equilibrium-state heterozygosity relative to 
its corresponding value for a single mate (s = 1). In con- 
clusion, multiple paternity promotes heterozygosity by pro- 
longing the duration of peaks of variation that reach the most 
distant islands. 



IV. DISCUSSION AND CONCLUSIONS 

In this study we analysed the effect of multiple paternity 
on genetic diversity over spatial and temporal scales in a 
metapopulation. We quantified the effect of multiple pater- 
nity during the colonisation of semi-isolated populations and 
in the equilibrium state developed after the colonisation phase. 
Our conclusions given below can be generalised to a metapop- 
ulation of patches that are partly isolated from each other, for 
example, by sandy beaches, harbors or other types of less suit- 
able habitats. 



We introduced a mating model which allows for different 
levels of multiple paternity in a population. In order to de- 
termine how realistic our mating model is, we compared the 
male family sizes of female broods obtained within our model, 
to empirically observed family sizes in populations of L. sax- 
atilis from natural habitats and under experimental conditions. 
We found that male family size distributions from the exper- 
iments, where the true number of mates was known, are in 
very good agreement with our mating model. The best-fitting 
parameters indicate fewer matings than the brood size, sug- 
gesting that some of the eggs are fertilised by sperm pack- 
ages retained between matings. By contrast, the correspond- 
ing empirical distribution in natural populations is best fitted 
by assuming an unlimited number of matings (i.e. no sperm 
retention). This is consistent with the high density of snails 
observed in the wild. However, the empirical distribution 
shows an excess of males with a single progeny compared to 
the mating model with the best-fitting parameters. This dis- 
crepancy could be explained by post-zygotic selection, where 
sperm from several males compete, resulting in uneven suc- 
cess among males. However, this effect should also be present 
in the experiments, but we did not find any evidence of it in 
the data. Another possibility is that this pattern is due to vari- 
ation in individual snails movements in the wild, where some 
snails might move around more extensively and mate with a 
new partner each time while others stay within a small area 
and mostly re-mate with individuals in the close surroundings. 
This variation in mating behavior cannot happen in the exper- 
iments where the snails are confined to each other. 

Within our mating model, by increasing the degree of mul- 
tiple paternity, the effective population size increases, and thus 
the heterozygosity increases. However, since mating is con- 
sidered to be costly H, it is not clear whether or not mating 
multiple males is an evolutionary strategy to increase the het- 
erozygosity. Recall that we estimated that in natural popula- 
tions of L. saxatilis the probability that two eggs are fertilised 
by the same sperm package is likely to be very small. Un- 
der our model, this probability is equal to zero if each sperm 
package fertilises one egg, or if the actual number of sperm 
packages a female receives during her reproductive cycle is 
very large. If the latter applies, we find that it is unlikely 
that the heterozygosity increase is the main reason for the 
extreme multiple paternity in this species. As earlier sug- 
gested, it seems likely that the cost of rejecting an intercourse 
is higher than the cost of accepting it, and a consequence of 
convenience polyandry Jgt]. Nevertheless, the consequences 
of multiple paternity for heterozygosity in relatively small and 
semi-isolated populations are substantial. This is summarised 
in the following. 

At a given distance from the mainland, populations with 
high degree of multiple paternity establish higher heterozy- 
gosity than populations with low degree of multiple paternity. 
While this effect is substantial during colonisation, it is mod- 
est in the equilibrium state. This is explained in the following. 
Upon the arrival of founder females to an empty island, the 
carrying capacity of an island is reached within a single gener- 
ation. Therefore, such a newly established population receives 
genetic material of most males that the founder females were 
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inseminated by. By contrast, in the equilibrium state, all moth- 
ers present in an island contribute to the population in the next 
generation, and hence the impact of immigrant females to the 
next generation is rather small. From this reasoning, we find 
that it is possible that the effect of multiple paternity upon het- 
erozygosity during colonisation might decrease if the growth 
rate of the island populations up to the carrying capacity were 
less than infinite (as assumed in our model). 

The heterozygosity at distances far from the mainland fluc- 
tuates significantly. Long periods of almost complete loss of 
genetic variation are interrupted by bursts of high heterozy- 
gosity, and this effect is boosted by multiple paternity. The 
durations of high- and low-heterozygosity phases could be an 
important survival factor in natural populations. For exam- 
ple, the low-heterozygosity phase could be disadvantageous if 
a malignant disease appears in the population, assuming that 
only a particular mutation (not present in the population, or 
being rare) can survive the disease. 

The wave-like nature of the spread of new alleles from the 
mainland population is also seen in the correlation of genetic 
diversity at neighbouring islands. We find that the correlation 
between heterozygosities at a pair of nearest-neighbouring is- 
lands increases as distance from the mainland increases (re- 
sults not shown). These results suggest intermittent bursts of 
genetic diversity in remote islands, an effect which becomes 
stronger as the degree of multiple paternity increases. 

The conclusions given above are confirmed by our com- 
puter simulations. In order to minimise computing time dur- 
ing simulations, we assumed that the mainland heterozygosity 
is equal to unity, which guarantees that whenever a migrant 
from the mainland comes to the first island, it carries genetic 
material that previously existed neither on the mainland nor on 
any of islands (and thus the population dynamics on the main- 
land does not need to be simulated explicitly). However, we 
emphasise that the conclusions given above qualitatively do 
not depend on the value of heterozygosity on the mainland. 

We note that, unlike in our model, it is possible that the 



rate of successful colonisation in natural habitats is smaller 
than the rate of migration. For example, if an immigrant fe- 
male carries a small number of progeny, her progeny alone 
may not be enough to colonise an empty island success- 
fully. By allowing for the rate of successful colonisation to be 
smaller than the rate of migration in the colonisation model, 
the equilibrium-state values of heterozygosity remain equal to 
those obtained under the assumption that the colonisation and 
migration rates are the same (as in our model). The values of 
heterozygosity during colonisation, by contrast, are expected 
to differ from those found here. 

In summary, this study can be used to quantify the gene 
flow between partly isolated natural populations using allelic 
frequencies at a number of neutral loci. Since our results show 
that the heterozygosity exhibits extreme fluctuations in popu- 
lations founded through repeated founder events, we raise the 
question of whether similar fluctuations can be observed at 
any given time at neutral loci sampled genome-wide. In order 
to answer this and related questions, the effect of recombi- 
nation needs to be analysed. Since island populations in our 
model experience at least one severe bottleneck, we expect 
that the degree of linkage disequilibrium in the colonisation 
phase is constant over a range of genetic distances, as shown 
in l20ll . However, how multiple paternity affects linkage dise- 
quilibrium during colonisation and in equilibrium is yet to be 
understood. It would also be interesting to analyse how selec- 
tion combined with recombination affects genetic diversity in 
a metapopulation. Such results would provide an advance in 
the endeavor of identifying genes under selection, and espe- 
cially, the genes underlying speciation l2lM23ll . 
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I. EFFECTIVE POPULATION SIZE UNDER MULTIPLE PATERNITY 



In this section we compute the effective population size under the mating model introduced in Section II of the main text. We 
assume that a population is diploid, isolated, well-mixed and that it consists of Nf females, and N m males. We assume Nf ^> 1 
and N m ^> 1. Note that in a diploid population with effective population size N c ^> 1, the population homozygosity F T in 
generation r is given by: 



F T &— e r + (l — )xr • (1) 

2N e V 2iV e V ; 

Here e T is the inbreeding coefficient which stands for the probability that two alleles within a single randomly chosen individual 
are identical at time r. The second term, % T , is the coancestry, that is the probability that two alleles sampled in generation r 
from two different randomly chosen individuals are identical. In what follows, we compute F T under our mating model, and 
thereafter we use Eq. (Q]i to determine the corresponding effective population size. We assume that mutation rate per generation 
per allele is fx <C 1, and we employ the infinite-alleles model []]]. Our calculations are based on the approach used in yj]. 

Under the assumption that all females are equally successful in producing offspring, it follows that the probability for two 
offspring to come from the same female, Pf, is equal to Pf = Nf . Let k be the probability that two offspring coming from the 
same mother also share a father. The probability that two offspring come from the same male, P m , is thus: 



P m = PfK+(l-P { )— , (2) 

where iV" 1 is the probability that two offspring having different mothers stem from the same father. The probability n has two 
contributions: the probability that two offspring come from the same sperm package, p, and the probability that two offspring 
do not come from the same sperm package but they come from the packages coming from the same male, (1 — p) / s. Here, we 
assume that sperm packages are chosen with replacement to fertilise eggs. It follows that k is given by 



P+{l-p)-, (3) 

s 



where p is 



UN \ 

P = TTlvM" ' (4) 

As mentioned in Section II of the main text, we take k _1 as a measure of the degree of multiple paternity. The case of single 
paternity corresponds to n = 1. 

Using the expressions for P m and Pf, we compute e T , and % T recursively. Under our model we find 
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e T = (1 - h) 2 Xt~i 
- ( - — 

+ 4\N~ f j 

After rearranging the terms in Eqs. (Ell-©, and by keeping only the leading-order terms, we obtain: 




(5) 



(6) 



e T = (1 - 2h)xt-i , 

1 fl + K 1 

Xt = ; 



1 fl + K 



1 



Ni N m J 8\ Nt 

2 fl + K 1 



1 



Xt-1 - 2/iXr-l 



8 V Mi 

Using Eqs. (|7])-([8]l we find the standard expression for the equilibrium-state homozygosity 



(7) 



(8) 



F = 



1 



1 + 6 ' 

where, as usual, 8 = 4fiN c is the scaled mutation rate. The effective population size A c is given by 



(9) 



A n =4 



1 + K 

N { 



1 

An 



(10) 



In the case N m = Nf = N (as assumed for island populations in our colonisation model), Eq. ( fT0T > becomes: 



N r , 



47V 
2 + k 



(11) 



Upon setting k = 0, Eq. (fTTT > reduces to A^c = 2 A, that is the effective population size becomes equal to the census population 
size. This is the maximum value that N c can acquire in an isolated population of census size 7V m + Nf, 



II. HETEROZYGOSITY IN THE COLONISATION PHASE 



In this section we compute the population heterozygosity in the colonisation phase. As mentioned in Section II in the main 
text, it is assumed that only the mainland is populated initially, and that its population size is constant in time. Moreover, we 
assume that the mainland population exists for long time before the start of colonisation (which occurs at generation r = 0). We 
denote the colonisation-phase heterozygosity on the mainland by Hc°\ 

(i) 

In order to find an expression for the heterozygosity of island i in the generation when it is colonised, H c 1 , we use the 
coalescent approach. Recall that a populated island is assumed to consist of N males and N females (N is large), and it is 
assumed that this population size is much smaller than that of the mainland. Moreover, in the following we assume that the 
migration rate M is small, M <§C 1. The average time between two successive founder events is thus long, and typically one 
founder female arrives at an empty island (the probability that two females come simultaneously is of order M 2 , and it is 
negligible for M <C 1). Under this assumption, the ancestral population size of the newly established population at island i can 
be represented by a sequence of i bottlenecks, such that each bottleneck lasts for one generation (since the founder female gives 
rise to 2 A offspring), and the time between two successive bottlenecks is on average generations long. Upon expressing 
the generation index r by t such that r = [2t N c \ , where A c is given by Eq. ( fTTT >. the waiting time between two successive 
founder events is approximately exponentially distributed with mean 
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{2MN c y 1 . (12) 

In order to compute the heterozygosity, we note that in our model, the mainland acts as the only source of genetic variation. This 
allows us to argue the following. First, if the MRCA of two alleles sampled randomly from the newly established population at 
island i was born on island j < i (j ^ 0), the two alleles sampled are identical. Second, if the MRCA was born on the mainland, 
the two alleles are expected to be identical with probability F^ ' = 1 — Hc°K Therefore, in order to compute Hi % \ it suffices to 
determine the probability that the MRCA of two lines sampled from the newly established population at island i stem from an 
allele that was born on the mainland, P(0\i). 

The probability P(0\i) has two contributions. The first contribution is the probability that the MRCA of two alleles sampled 
at island i is not found during a bottleneck. We find this to be equal to 1 — |(1 + «)■ The second contribution is the probability 
that the MRCA of two alleles is not found between two successive bottlenecks. This term is equal to 2MN C {2MN C + 1) . It 
follows that P(0\i) is given by 



8' ') \2MN C + 1 
Therefore, the colonisation-phase heterozygosity at island i is: 



iT« = P(0\i)H^ . (14) 

Note that for the case described here, the population size switches between 2N (N males and N females during the waiting time 
before the colonisation of the next island) and unity (one inseminated mother during a bottleneck). Therefore, the probability k 
does not enter Eq. (fl4l only through the expression for N c , Eq. ( fTTT i. 

The heterozygosity in the colonisation phase for different parameters of our model is shown in Fig. Q~K-C. The analytical 
result, Eq. (fT4T >. is shown as solid lines, and the results of our computer simulations are shown as symbols (the solid lines in 
Fig- 03* correspond to the solid lines in Fig. 4A in the main text). We see that the agreement between Eq. ( TBI and the results of 
computer simulations is good for M = 0.05, whereas for M = 0.5, Eq. ( fT~4T > underestimates the results of computer simulations. 
This is discussed in Section III in the main text. 



III. HETEROZYGOSITY IN THE EQUILIBRIUM STATE 

In this section, the equilibrium state within the model introduced in Section II in the main text is analysed. Expressions for 
the equilibrium-state heterozygosity at distance i = 1, . . . , k from the mainland are derived under the assumption that all islands 
are populated (that is, the colonisation phase is over). As discussed in SectionQ]of this supplementary material, the inbreeding 
coefficient in generation r at island i, e£ l , and the coancestry Xt contribute to the homozygosity Ff at island i in generation 
r. The coancestry between islands i, and j is equal to the inter-island homozygosity, and for this case we use the notation 

As mentioned in Section II of this supplementary material, we assume that the mainland is the only source of genetic variation. 
All habitats are assumed to have equal numbers of males and females. The population size on the islands i = 1, . . . , k is assumed 
to be large (and equal to 2N), but much smaller than that of the mainland. As before, the heterozygosity on the mainland in 
generation r = is denoted by ■ 

According to the spatial model introduced in Section II of the main text, the female-migration rate per island per generation 
is 2M for islands i = 1, . . . , k — 1, whereas for the mainland and for the island furthest from the mainland it is equal to M. 
Since the population size on the mainland is much larger than that of a populated island, the process of migration does not affect 
genetic variation on the mainland. Therefore, we have H T ' = i?c°' ■ For the island populations, we consider separately the 
inter-island and the intra-island homozygosity. First, we treat sampling from two distinct islands, i ^ j. Second, we consider 
the case of sampling within a single island i = 1, . . . , k. Our calculations given below are based on the approach employed in 

i. 

The inter-island homozygosity i^Yy for i = 0, < j < k satisfies the following recursion: 

Ft°+i = (1 - m + S^Fpfi + | (F^-V + (1 - S^FW+Q) . (15) 
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FIG. 1: Heterozygosity in the colonisation phase (A-C), and in the equilibrium state (D-F). The results corresponding to Eq. d 1 4b are shown as 
solid lines, and the results of computer simulations are shown as symbols. The number of available mates is s = 1 (blue), s — 2 (red), s — 3 
(green), s = 5 (magenta), and s — 10 (black). Averages are over r = 10 4 independent realisations of the process of colonisation of empty 
islands in A, that is over r — 2 ■ 10 4 realizations in B, and C. In panel D, averaging is done over r = 1.5 ■ 10 7 generations (the initial r = 10 7 
generations being discarded), and over three independent realisations of our model. In panel E, averages are over r = 4 ■ 10 7 generations (the 
initial r = 5.5 • 10 7 generations being discarded) and over four independent realisations of our model. In panel F, averaging is done over 
r = 5 • 10 7 generations (the initial r = 5 • 10 e generations being discarded) and over five independent realisations of our model. Remaining 
parameters used: mainland heterozygosity = 1, number of islands k = 10. 

Here m = 2M/N <C 1 is the migration rate per island per female per generation, and Sj k is equal to unity when j = k, and it 
is zero otherwise. The inter-island homozygosity for < i < k, < j < k, i ^ j, obeys: 



(l- m )jrW)+ mx W fe. i+^-.i) 



+ |(i-^._ 1 )(f^'- 1 )+f^)) 



+ 0(m 2 ) 



Lastly, when i = k, < j < k, we find 



(16) 



(k,j) 



T+l 



(l-^)(l- m) F(M 



+ y (1 - y)(^ +1) (l - S W ) + 4 J+1 xt fej+1) ) 

+ ^(i- m )(Ff-^(i-4 J+ i)+xr ij) ^ + i)+o( m 2 ). 



(17) 
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The inbreeding coefficient of the population on island < i < k satisfies: 

e«=(l-m) X « + ^- 1) + ?xr i) , (18) 



and the coancestry is given by: 



XXL = (1-mV— 

At+1 v ' N(l-m) 



+ m(l - to) (i^ 1 ) + F^- 1 - 1 )) + 0(to 2 ) . (19) 
For the island furthest from the mainland (i = k), we have: 

a, = f i - x< fe) + ^ , (20) 



(jo = (1 _ ™ )2 i / l+iz^n + K) + + x^!\ 

Xr+i U 2 ) N{1 _™ ) y 8 U + M + 4 X T + 2 j 
+ (i - ™) 2 (i ( V + + (1 - 1)^ 

+ m(l-y)F^- 1 )+0(TO 2 ) . (21) 

In what follows we keep only the leading order terms in Eqs. ([T5b-(|2T]). Moreover, we use the scaled time t, where r = [2 A?" e tJ , 
and N c is the effective population size of an island population under our mating model (see Eq. (fTTI) ). In these units of time, we 
denote the homozygosity at time t at island i by (t). 

We finish this section by giving differential equations (with only the leading order terms) for the mainland-island homozygos- 
ity, then for the inter-island homozygosity and, lastly, for the intra-island homozygosity. 

For the mainland-island homozygosity we find: 

= -dtF (0 ' l) {t) + M c ^F {0 - i+1} (t) + F (0 ' l - 1} {t)-2F {0 - l} {t)^ . (22) 
Here, M e = 2M N c /N, and it is assumed that i < k. For i = k, we have: 

= -d t F {0 ' k \t) + Me (F(°' k -V(t) - F(°' fc >(t)) . (23) 
For the inter-island homozygosity between islands i and j, where < i < k, < j < k, j ^ i, we find: 

= - d t F^\t) + M c (1 - 5^) (F^+^it) + F^-^it)) 

+ MA-u (F^it) + F^i)) + M c (1 - 5i+ ld ) (F<M-V(t) + F^ 1 -^)) 
+ M c 5 l+l ,j (f W (<) + F (i+1) (<)) - 4M F (i J) (t) . (24) 
For i = k, < j < k, we obtain: 
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= -d t F {k ^ (t) + M e (t) - 2T (fej) (<)) . (25) 

Finally, for the homozygosity at distance < i < k from the mainland we find: 

= (-ft - 1) F« (t) + 2M e (i^ 1 -*) (t) + F^- 1 '*) (t) - F« (t)) + 1 , (26) 
For the island furthest from the mainland, i = k, the corresponding expression is: 

= {-d t - 1) (t) + Mo (F (fc ' fc - 1} (t) - 2F (fc) (t)) + 1 . (27) 

By setting in Eqs. d22b-(|27l) the terms involving the time derivative to zero, one finds the expressions for the equilibrium- 
state homozygosity of the system. The equilibrium-state heterozygosity is obtained upon subtracting the equilibrium-state 
homozygosity from unity. Upon setting H^°> = 1, the results shown in Fig.QJ)-F are obtained. We note that the lines in Fig.QJi 
correspond to the lines shown in Fig. 4B in the main text. 

IV. DURATIONS OF LOW- AND HIGH-HETEROZYGOSITY PHASES 

A scheme for computing the durations of low- and high-heterozygosity phases is shown in Fig. [2}V. As indicated in this 
figure, we consider values of the heterozygosity smaller than 0.1 to be low, and values of the heterozygosity larger than 0.4 to 
be high. The maximum value for the low phase (0.1) is chosen because a locus is commonly considered monomorphic if the 
heterozygosity at this locus is < 0.1 (see 14J1). The minimum value for the high phase (0.4) is chosen since the typical maximum 
value that the heterozygosity has at the island furthest from the mainland is ss 0.5 for the parameters chosen in Fig. [2}V. Note 
that the maximum value of the heterozygosity at a locus with only 2 allelic types is equal to 0.5. 

The method used to record the times of transitions between low- and high-heterozygosity phases in computer simulations is 
explained next. Say the heterozygosity is in the high phase at time r = 0. We record first the nearest point in time when the 
heterozygosity becomes less than 0.1. Say this occurs in generation To. Second, we search for the last generation before To in 
which the heterozygosity has a value larger than or equal to 0.25 (the middle point between 0.1 and 0.4). Say this happens in 
generation t\ < tq. We take t\ to be the time of transition from the high- to the low-heterozygosity phase. The transitions 
from the low to the high phase are recorded using a similar scheme. The durations of the high- and low-heterozygosity phase, 
Thigh and Ti ow , relative to their corresponding values for s = 1, computed as explained above, are shown as symbols in Fig.|2jJ, 
D, and F. For comparison, we show in Fig. [2]C, E, and G the corresponding equilibrium-state heterozygosities relative to their 
values for s = 1 . 

Next, we briefly explain how to estimate 21 ow , and Thigh analytically. The heterozygosity can switch from the low to the high 
phase if new genetic material comes to the population, and if this material is not lost due to random genetic drift. Recall that, 
in our model, migration is the only process bringing new genetic material to islands. However, some migrations bring genetic 
material that already exists in a given island population (but note that when H^°> = 1, then the first island receives new genetic 
material with each migration). Yet, it is possible to estimate the effective successful migration rate per allele per generation, 

(i) 

Too , at a given distance i from the mainland using the analytical results derived for the equilibrium-state heterozygosity at 
island i, namely ffW. Here, the term 'successful' implies that migrants bring new genetic material to the population. Note that 

(i) 

under the assumption that a new allelic type is introduced to the population at distance i from the mainland at rate to c per allele 
per generation, the equilibrium-state heterozygosity is computed as 



H«. . , ,28) 

1 + 4mi' ) JV c 

(i) 

where 2to c N c is the total number of new allelic types introduced in the population per generation (this expression is analogous 
to the usual expression for the heterozygosity which involves the scaled mutation rate, see Eq. ©). In the case m^' N e <C 1, 
it follows that typically every (2mi^ A^) -1 generations, one mother comes to an island carrying a new allele. Furthermore, it 
is known that a large haploid population with effective size 2A r c spends on average 2/N c generations in the state with this 
allele having the frequency N c before fixation occurs. It follows that the typical waiting time for such a successful establishment 
of new genetic material at island i is equal to (4m e This estimate agrees well with the results of our computer simulations 
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(results not shown). 

We now explain how to estimate the average duration of the high-heterozygosity phase, Thigh. in the case of rare income of 

(i) 

new allelic types to island i, rrie N c -C 1, the island population at a given locus typically has at most 2 alleles. Thus, Thigh can 
be approximated by the time that a locus with two alleles in a Wright-Fisher population of N c diploid individuals needs to reach 

(i) 

fixation. Note that under the condition m e N e <C 1, fixation occurs typically much before new genetic material arrives to the 
population. The number of generations until fixation in a diploid population at a locus with two alleles, ti oss , is J3l 

ri oss (a) ss -4N e [a ln(a) + (1 - a) m(l - a)} . (29) 

Here a is the initial frequency of a given allelic type. For a given value of the number of available mates, s, we compute Thigh 
upon integrating Eq. ( 1291 from a = 0.1 to a = 0.9. The integral boundaries correspond to the value of heterozygosity w 0.2, 
which is close to the value of 0.25 in the method used for determining the transition time between the low and the high phase. 
We note that Eq. d29l results in underestimated time of fixation if the rate of income of new genetic material is not too small. 
This is confirmed by the results shown in Fig.|2jj, D, and F. Note that the solid lines in Fig.|2)i, and C correspond to the lines in 
Fig. 5B, and C (respectively) in the main text. 
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FIG. 2: (A) Illustration of the method used to determine the duration of low- and high-heterozygosity phases. The heterozygosity represented 
in terms of the low and high phases is shown by the magenta line. The black line depicts the result of computer simulation. The data shown 
correspond to those in Fig. 5A in the main text, for island 10. (B), (D), and (F) Durations of low- and high heterozygosity phases relative to 
their corresponding values for s = 1 (blue and red, respectively). (C), (E), and (G) Equilibrium-state heterozygosity at island 10 relative to its 
corresponding value for s = 1. The results of computer simulations are shown as symbols, and the analytical results are shown as solid lines. 
The parameters in B, and C correspond to those in Fig.QJ). The parameters in D, and E correspond to those in Fig.QJi. The parameters in F, 
and G correspond to those in Fig.[]J\ 
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